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SUMMARY 

Operators  of  numerical  integration  of  functions  are  constructed. 

Also  considered  is  the  question  of  estimating  the  total  error.  The  operators 
are  applied  to  the  solution  of  a  problem  with  initial  conditions  for  ordinary 
linear  differential  equations  with  constant  coefficients.  It  was  determined 
that  the  process  of  finding  the  solution  of  -a  differential  equation  of  the  type 
being  examined  by  means  of  numerical  integration  operators  reduces  to  arith¬ 
metic  operations  with  numerical  matrices  and,  consequently,  is  highly  suitable 
for  computer  realization.  An  example  is  given,  illustrating  the  advantages 
of  the  proposed  method. 
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19.  ABSTRACT 


Operators  of  numerical  integration  of  functions  are  con¬ 
structed.  Also  considered  is  the  question  of  estimating  the  total  error.  The 
operators  are  applied  to  the  solution  of  a  problem  with  initial  conditions  for  ordinary 
linear  differential  equations  with  constant  coefficients.  It  was  determined  that  the 
process  of  finding  the  solution  of  a  differential  equation  of  the  type  being  examined 
by  means  of  numerical  integration  operators  reduces  to  arithmetic  operations  with 
numerical  matrices  and,  consequently,  is  highly  suitable  for  computer  realization. 
An  example  is  given,  illustrating  the  advantages  of  the  proposed  method. 
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3Y  MEANS  OF  NUMERICAL  INTEGRATION  OPERATORS 


A.  V.  Nesterchuk  [Nestercuk] 


Considered  in  [l,2j  are  several  aspects  of  the  solution  of  differential 


equations  by  means  of  numerical  integration  operators.  Considered  in  the  present 


paper  are  new  numerical  integration  operators  and  their  application  to  the  computer 


solution  of  ordinary  linear  differential  equations.  Also  examined  is  the  problem 


of  estimating  the  total  error. 


Given  a  function  f(x)  continuous  on  the  segment  [0,lj.  We  divide  the 

segment  [ 0, l]  into  n  equal  parts,  h  «=  jk  We  introduce  the  notation 

x0  —  0,  x,  ~  ih,  ft  —  f  (tj),  |/jl  =  l/o*  •  >  •  «  f„b  /  ~  0, 1 , 2 . tu 


We  shall  calcuiateJ7(/)df,  0*x<lfrom  one  of  the  quadrature  formulas, 
e.g.,  by  the  trapezoidal  formula.  Then 


* 

0 

J /  (Qd/ssA-JLtf0  +  2f,  +  fa)  =  h  /,)•  -  p  j 


Translated  from  Ukr.  Hat.  Zh.  [Ukr.  Math.  J.3»  Vol,  17,  No.  4,  pp.  112-119 
<1965). 


(1) 


n 


A* 


f  (t)  dt  %  ft. 1  tf0+  2/f  4  . . .  +  2f„_t  +  /„)  =  ft  l/e.  /, . Ui 


2 

1 


or,  in  the  general  case. 


where 


iHOdt  «| fj.fc.il,,  ft  =  0.1, 2., 

u 


•  )i|  i  ~  0^  ly  2,  •  *  •  t  /I* 
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r  o  i  i  i .  . 

0122.. 

0  0  12.. 

0  0  0  1.. 
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2 

2 
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0  0  0  0  ...  1 


Thus,  the  matrix  is  a  difference  operator  on  the  right  in  the  trape¬ 


zoidal  formula.  In  what  follows  we  shall  call  it  the  numerical  integration  oper- 

th 

ator.  Denoting  jf(i)dt,  equality  (1)  can  be  rewritten  in  the  form 


ft  =  0,1.2,..../;  £  =  0, 1,2 . n. 


Here  the  error  can  be  found  on  the  basis  of  the  equalities 

=  +  ‘  0  <*!<*• 


y3  =  +  2f,  +  fj  -  4'  W  (5,),  0  <5,  -4  2ft, 


(2) 


+'*»  +  — +#_  +W— 0 


From  (2)  it  follows  that 


(3) 


where 


ftfft* 


=  -  j2*  [0, 1. 1. ... .  IJ,  M  «  nwx  \nb  |. 


o<ur 


The  one-row  matrix  will  be  called  the  error  vector.  Expression  (3) 


characterizes  the  error  of  the  method.  The  total  error  is  the  sum  of  the  method 
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error  and  the  computational  error.  Apfopos  the  total  error,  the  following  theorem 


holds. 


Theorem  1_.  I_f  the  f unct ion  f (x) ,  defined  and  continuous  on  the  segment 
LO.l],  has.  on  this  segment  a  second  derivative  for  which  the  following  condition 


is  fulfilled: 


A 

then  in  calculating  J  0  <  1,  b^r  means  of  the  operator  A^,  when  the  errors 

o 

in  the  calculated  values  of  the  function  f*(ih)  (i  =  0,1,.. .,n)  do  not  exceed  the 

jjj.2 

quantity  ,  the  vector  of  the  total  errors  R*  has  the  form 

10,2,3 . *+M. 

Proof.  Let  f(x)  be  the  exact  and  f*(x)  the  calculated  values  of  the 
function.  The  total -error  vector  for  the  approximate  vector  h[f*]A^  is  deter¬ 


mined  in  the  form  of  the  difference 


JJttOdf  — MPM. 


AtsO,  1, 


We  calculate  the  i-th  component  of  vector  (4).  Taking  (2)  into  account,  we  get 


~/t(f0 +2 2f,_, + f ,) -j2htrd)~h\rvfy . r,\~j 


Yh  tfo  +  2/,  +  •  •  •  +  2  -  1ft  (f;  +  2r,  + . . . 

•  •  •  +  2/J_,  +  /J)  =  ^rh  {iff—  fp  -j-  (2 f,  —  2fJ)  -f . . . 


•  •  •  +  2  (f._, + (f'-ffl—Lhr®, 

i  =  0,  l,2,...,n;  0<|  </ft. 


By  the  condition  of  the  theorem  | f  -  fj  £  yy-  »  then,  from  (5) 


3- 


(6) 


( I'w<"  ~*l04,)J<T1*Tr,+n 


mt 

—  ~I2~{£  +  1).  i«=»  1,2,3, . n. 


Taking  the  obvious 


equality  J([J /W*J j  _0  and  (6)  into 


account,  we  can  write 


—  ~|2~  1®*  2*  3,  • .  • ,  n  *}*  1 1- 


More  accurate  will  be  an  operator  constructed  by  means  of  second-  and 
third  degree  parabolas.  When  x  *  h  we  use  the  trapezoidal  formula;  then 

H 

V 

When  x  *  2h,  by  the  Simpson  formula  we  have 


34 

Jf  (0*^- /«(/„  + 4/, +;> 


and  when  x  **  3h,  we  have  by  the  "3/8"  Simpson  formula 


» 

J  f{f)dt  sgl -fc(f#+  3/t  +  3/,  +  fj). 


But  if  x  *  4h,  we  apply  the  "one-third"  Simpson  rule  twice,  acid  when 
x  »  5h  tnis  same  rule  is  applied  from  0  to  2h  and  the  "three-eighths"  Simpson 
rule  from  2h  to  5h.  As  a  result,  by  analogy  with  (1)  we  get 


h[ft I<4j»  4«0,l,2,...,i;  *  =  0, 1, 2, , 


where 


0,-5^  r;j^.t«. f  ?.;* ': F. -:,' >> ' ^~'"-.-tJ *' ?  ' -  I'^’Ti -^xl- ?»i£~S'i ^~}C ■xS’^'c  I'' .—  '-  ?  ’  _'”Ss‘,9!?y5s*4-,-i'^, 


ro  12  8  s  8  8.  ..8  - 

0  12  32  27  32  32  .  .  .  32 


0  0 
0  0 
0  0 
0  0 


8  27  16  17  ...  16 

0  9  32  2?  ...  32 

0  0  8  27  ...  16 

0  0  0  9  ...  32 


LO  0  0  0  0  O  ...  8(9) 


The  last  diagonal  element  equals  eight  if  n  is  even  and  9  if  n  is  odd.  The 
error  in  this  ease  can  be  found  on  the  basis  of  the  equalities 

0,  ~T*(f*  +  f,)-j2mi,h  0  <it  <h, 

0i®|-Atf9  +  </,+/^*~^A*/(IV,{|}).  0<£a<2A, 

ys  -  ~  A  (/, 4-  3f,  4-  3/j  +  /,)  - 1  hSf™  «,),  «  <  *,  <  3ft. 

£4=±A(/e4.  4/>  +2/i  +  4?* + '<) +  fmxG$' 

0<%<2h,  2h<$<4ht  (] 

06  =  *<8/e  +  32/,  +  >7/,  +  27/,  +  27/,  +  9/j)  - 

-^•(3/^^  +Zlfav'(Q)> 

0<|;<28,  2A<K<SA, 


Thus,  it  follows  from  (11)  that 

nw-AifMi***  <l2> 

where 

«,  =.—•  (0,  t.  8. 27. 16. 35. 24. . . .], 

Af  =  max(i60ft-*r(£){t  f/fW»(6)|fc 

o*$d 

ft  =  0,  l,2,...,i:  i=sO,  I,2,.,v,n. 


Insofar  as  the  total  error  is  concerned,  a  theorem  analogous  to  theorem  1 


is  fulfilled  in  the  process  of  applying  the  operator 

Theorem  2.  If  a  function  £ (x) ,  given  and  continuous  on  the  segment 
L 0, l] ,  has  on  this  segment  derivatives  of  up  to  and  including  the  fourth  order 
£Sil£  the  conditio-.!  |6?/j-'I/*{$)|  < M,  {f'Kt(l)\<M.  0  <£<  1,  are  fulfilled,  in 


r 

calculating  J/(f)<£»  bj?  means  o_f  the  operator  when  the  errors  of  the  cal¬ 


culate  1  values  of  the  function  f*(x)  in  each  integration  step  h  do  not^  exceed  the 


v»-lue  ,  the  total-error  vector  R."  has  the  form 


^  "  m (0* 2* I0<  20, 40, 30, 50, , . 


lha  <  oof  of  theorem  2  Is  analogous  to  the  proof  of  theorem  1,  and  there¬ 


for"  ve  viil  r  /e  it  here. 


Up'-.-  tepeated  integration  by  means  of  numerical  operators,  the  matrices 


representing  these  operators  multiply  out.  In  fact,  if  ve  calculate  ,?  J , 

\dt\f(x)dx, 


I 

then,  denoting  J/(T)rfr  #(t),  we  get  from  (1) 


««  1,2, 


It  is  easy  to  see  that 


jdfj  *♦••£/(«)<*“  ft" 


t  =  i,2. 


It  should  be  noted,  however,  that  it  is  not  always  expedient  to  use 


formula  (14),  since  the  calculations  can  be  made  more  successfully  if  the  left- 


hand  part  of  formula  (14)  is  replaced  by  the  Cauchy  formula: 


w  n,  »i« 


fdn,  j\f%  l dnr..  j*  =  Ci5) 


and  then  formula  (1)  is  applied. 


Now  let  us  consider  the  solution  of  ordinary  linear  differential  equa¬ 


tions  with  constant  coefficients  by  means  of  numerical  integration  operators. 


£k 


W 
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Given  the  differential  equation 


Sid 


(16) 


(17) 


Qy' •rby**f{x),  Q  <x  <  1, 

^ere  a  and  o  are  constants,  with  the  initial  condition 

t/{0)  =  y0. 

By  numerical  integration  of  Eq.  (16)  taking  (I)  and  (1?)  into  account  we  gee 

a*;,  ji.  5 . 5}  •\-bh\yk\Al-  h]f)}Al 

or 

{£*1(3/ +  &H)  ~h\fJAs  +  ay0l't,i . 15.  »«1. 2.  (18) 


where  I  is  the  unit  matrix.  The  raatrix  al  +  bhAj^  is  invertible $  since  it  is 


triangular  and  its  diagonal  elements  are  different  from  zero;  and  when  h  ?  -  r-~- 

b3jj 


(a.,  are  the  diagonal  elements  of  the  raatrix  Aj)  det  (al  bhA^j  r  0.  From  (18) 


l'yaI  =  mt]At  4-  ay0ll,  I, . . . ,  1JJ  B~\ 


(19) 


where  B  =  al  +  bhA^  (i  *  1,2,),  is  a  matrix  which  we  will  call  "resolving”. 


Equality  (19)  is  the  computational  formula  for  finding  the  solution  to  the  problem 
(16)- (17). 

Let  us  now  consider  the  second-order  differential  equation 

«/  +  by'  +  cy  =  /(*),  x*  |0, 1J.  (20) 


y(0)  =  ^  y1  (0)  =  $rj, 


where  a,  b,  c  are  constants. 

Integrating  the  first  time  we  get 


a{y‘  —  y$  +  yj  4- cj  ydt  =zjf(/)d/. 


(21) 


*+*+'1*-  r**+«+* 

Now  we  integrate  (21)  taking  (13)  into  account  and  get 

oy  +  b$ydi  +  c  [dlj  ydx  =  j  dl  |  f  (f )  rft  +  J  (ayt  4-  6ue)  di  -f  ay0, 

-7- 


2*  -K-^* 


or 


(22) 


iy*l  +  th  ff/*l  4,  +  CA  •  [^1  y4*  r,  h2  [/,]  A)  + 

+  Wo  +  At/^10, 1,2 . n)  +  fly0{l,  1 . 1),  t  =  i#  2; 

k  =  0,  l,2,...,n. 

For  the  given  equation  we  call  the  matrix 

B  — a!  +  bhAi  +  ch2A2,  i  =  1,  2. 

the  "resolving"  matrix. 

In  view  of  the  fact  that  the  matrix  is  triangular,  for  its  invertibil- 
ity  it  is  necessary  to  choose  h  such  that 

a  -}-  bhan  +  ch2dln  «?fc  0, 

then  det  B  r  0  and  formula  (22)  can  be  written  in  the  form 


fyj  =  {^lfk\A}+  (ay‘a  +  by0)  JO,  1 , 2 . n]  +  ay^X,  1 . 1J}  B~\ 

*=•'!,  2;  k  =  0, 1 , 9. . n. 


(23) 


Formula  (23)  will  be  the  computational  formula  for  finding  the  solution  to 
problem  (20).  When  integrating  Eq.  (21)  it  is  also  possible  to  use  equality  (15) 
instead  of  (13).  We  then  get 

(#*!  [at  +  bhA,  +  ch7  Jfc  —  /)„  = 

=  h2  [/tl  [k  —  j]iAi  -f-  {au0  +  by0)  (0, ! ,  2, . . . ,  n\  +  ay0  [1 , 1 . 1 ), 

where  i  -  1,2;  j  «  0,1,2,. ..,k:  k  =  0,1,2, ...,n;  [k  -  j]d  is  a  diagonal  matrix. 

We  proceed  analogously  in  the  solution  of  an  m-order  linear  differen¬ 
tial  equation,  applying  m-fold  integration  using  equality  (14)  or  replacing  by 
formula  (15). 

From  (19)  and  (23)  it  is  clear  that  the  process  of  finding  the  solution 
of  a  differential  equation  reduces  to  arithmetic  operations  with  numerical  matrices 
which  makes  this  method  suitable  for  use  in  the  computer  solution  of  differential 
equations  of  the  indicated  type,  since  standard  programs  are  used  to  multiply 
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and  invert  matrices. 


Here  the  process  of  inverting  the  matrix  B  has  an  appreciable  effect 

on  the  accuracy  of  calculation.  This  must  be  done  by  Che  method  indicated  in  [3], 

*  | 

i.e.,  by  introducing  a  correction  to  the  elements  of  the  matrix  B  . 

Example.  Let  us  consider  the  simplest  differential  equation 

=  —  4*.  x  e  JO,  1]  (24) 

with  initial  condition  y(0)  =  1,  having  the  obvious  solution  y  =  4x^  - 
6x  +  3  -  2e~^K.  For  the  solution  we  use  the  operator  with  h  =  0.2. 
The  computational  formula  for  finding  the  solution  of  Eq.(24)  has  the 
form 

(0.2^1^+ [i.i . ljjiH, 

where  £(x)  =  8x2  -  4x,  B  =  I  +  0.4A,,,  k  =  0,1, 2,3, 4, 5. 

The  solution  matrix,  obtained  analytically,  is  written  as 
[ykJanai  =  [1,  0.6194,  0.3410,  0.2376,  0.3562,  0.8343] ; 
obtained  by  the  proposed  method  it  is  written 

[y.  ]A  «£i  0.6205,  0.3408,  0.2365,  0.3552,  0.8331]; 

K  &2 

and  by  the  Runge-Kutta  method  (with  computational  formula 

Ay,  w  4  v  <*,  +  *j).  =  ft/  <*i«  ft)-  kt  “*/ (*/  +  ?*•  ft+y  A,) .  *#  =  V  (*,■  +  h>  ft  +  2ta  ~  *i)> 

[y,1  "IX  0.6139,  0.3369,  0.2334,  0.3512,  0,7246]. 

K  RK 

The  relative  errors  at  selected  nodes  are,  respectively: 
by  the  proposed  method: 

bl  «  0.167.,  52  *  0.077.,  &3  -  0.46%,  5^  =  0.287.,  i>5  »  0.14%; 
by  the  Runge-Kutta  method: 

TL  «  6.897.,  »  0.677.,  T),  =  1.77%,  tl  -  1,57%,  0  *  13.2%. 


Setting  %  -  z  in  (13),  we  get 

~  V  (Q/  (2)  +  P,(2)M  (z))*  [Q,  (2)  +  P,  (2)  M  (2)1 . 


2  —  2 


(mO 


(14) 


Starting  with  relations  (9)  and  (14),  and  then  completely  repeating  the  proof  of 
lemma  4.1  from  [l],  we  arrive  at  the  following  statement. 

Lemma  3> .  In  the  completely  indefinite  case,  for  every  finite  domain  G 
of  the  complex  plane  there  exists  a  constant  such  that 

Jq^Q^z)  <qQ.  (15) 

/-a 

Let  us  give  without  proof  one  more  statement  which  describes  the  M(z) 
at  any  fixed  point  z  (Im  z  /  0)  that  correspond  to  all  orthogonal  operator  distri¬ 
bution  functions  t(X)  of  the  positive-definite  function  F(t). 

Theorem  In  the  completely  indefinite  case,  for  fixed  z  (Im  z  f  0) 
the  set  of  M(z) 


A*(2)~0  (2)  +  j-U  (2)  Ur'*  (2) 


(16) 


corresponding  to  all  possible  orthogonal  distribution  functions  " (>.)  of  the  function 
F(t)  ranges  over  the  "operator  circle11  Kz  [l]  with  center 


wm  m 

i£+^PJ(2)Q,(2) 

/•o'  /» 0  . 


and  radii  (left  and  right)  rn  «=  - — L=r  r,/a (5) ,  r,(  =*  j — ~  rl/3(z),  U  ranges  over  all 

unitary  operators  in  H. 

Replacing  z  by  %  and  %  by  t  in  (13)  and  multiplying  we  get 
M  (t)  £,  ft,  2)  +  £,<$,  2)  =  |M  (2)  De  (6,  z)  +  Dt  (l,  2)J  M  {$), 


(17) 


where 


2) — p;(i)P/(E).  £>, «,  2>=£-a-2)^  qjwp,® 


/-6 


£0 *)-*-+ e  -  2)  y  p;  w  q,  ft).  £,(*.*>-  «  -*)  £  %<*)  q,  ft)  ■ 

/«o  /• 0 


(18) 
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